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ABSTRACT 

We present new results from a series of ACDM simulations of cluster mass halos re- 
solved with high force and mass resolution. These results are compared with recently 
published simulations from groups using various codes including PKDGRAV, ART, 
TPM, GRAPE and GADGET. Careful resolution tests show that with 25 milhon 
particles within the high resolution region we can resolve to about 0.3% of the virial 
radius and that convergence in radius is proportional to the mean interparticle sepa- 
ration. The density profiles of 26 high resolution clusters obtained with the different 
codes and from different initial conditions agree very well. The average logarithmic 
slope at one percent of the virial radius is 7 = 1.26 with a scatter of ±0.17. Over the 
entire resolved regions the density profiles are well fitted by a smooth function that 
asymptotes to a central cusp p oc r~'^ , where we find 7 = 1.16 ± 0.14 from the mean 
of the fits to our six highest resolution clusters. 

Key words: methods: N-body simulations - methods: numerical - dark matter — 
galaxies: haloes — galaxies: clusters: general 



1 INTRODUCTION 

A highly motivated and well defined problem in 
computational astrophysics is to compute the non- 
linear structure of dark matter halos. This is espe- 
cially timely given the abundance of new high reso- 
lution d ata that probe the central structure of galajc- 
ies (e. g. !de Blok et al." 20018; do Blok. McGaugh fc Rubir 



2001bl: FMcGaugh. Rubin fc do Blok 200 L: Swaters et al 



20031: Ide Blok fc"Bosmall2002l: iGentile et al]l2004D and clus- 
ters (e.g. ISand et al J 12004) . Furthermore, a standard cos- 
' mological paradigm has been defined that gives a well de- 
fined framework within which to pe rform numerical ca l- 
culations of structure formation (e.g. ISoergel et al.l 120031) . 
This subject has developed rapidly over the past few years, 
building u pon the pioneering results obtaine d in the early 
1990's bv iDubinski fc Carlbcrg (1993) and IWarren et. all 
il992h . More recently, the systematic study of many ha- 
los at a low resolution led to the proposal that the profile 
of an 'average' cold dark matter halo in dynamical equi- 
libriu m could be fi t by an universal two parameter func- 
tion jNavarro. Frenk fc White 1996), with a slope of that 
asymptotically approaches — 1 as r 0. At the same 
time, the stu dy of a few halos at high resolution questioned 
these results llFukushige fc Makindll997l: iMoore et alJll998t 
iMoore et alJll999t Ijing fc Sutoll2000t iGhigna et al.ll200Cl) . 
These latter authors claimed that of the order a million par- 
ticles within the virialised region where necessary to resolve 
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the halo structure to 1% and the slopes at that radius could 
be significantly steeper. Just within the last few months, we 
have seen several groups publish reasonably large samples 
of halos simulated with the necessary resolution that we can 
finally determine the sc atter in the density profiles across 
a range of mass scale s fFukushi ge. Kawai fc Makino 2004 

yw 



Tasitsiomi et alJl2004l: J^^nbsganss^B^^^^Os^ker^j2004 



Havashi et al.ll2003l:lNavarro et alJl2003HReed et. aJl2003l) . 

Much of the recent controversy in the literature has 
been due to limited statistics and the lack of agreement over 
what is a reliable radius to trust a given simulation with 
a given set of param eters. Several stud i es have attempted 
to address this issue llMoore et a i]|l998t iKnebe et al ]|200d : 
iKlvpin et al]|200ll : IPower et al.ll2003l:lDiemand et alJl20oll . 
Integration and force accuracy can be understood using con- 
trolled test simulations. However, discreteness is probably 
the most important and least understood numerical effect 
that can influence our numerical results which is exacer- 
bated due to the lack of an analytic solution with which to 
compare simulations. Our particle sampling of the nearly 
coUisionless fluid we attempt to simulate can lead to energy 
transfer and mass redistribution, particularly in the central 
regions that we are often most interested in. 

CoUisional effects in the final object or in the early hier- 
archy of objects can be redu ced by incr easing the number of 
particles A' in a simulation jPiemand e t al. 2004|) . The lim- 
itation to the phase space densities that can be resolved due 
to discreteness in the initial c onditions can also be overcome 
by increasing the resolution jBinnevll200l . As we increase 
the resolution within a particular non-linear structure, we 
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find that the global properties of the resolved structure is 
retained, including shape, density profile, substructure mass 
functions and even the positions of the infalling substruc- 
tures. This gives us confidence tha t our A'^-body calculation s 
are not biased by using finite N jBaertschiger et al.ll20o3) . 
The fact that increasing the resolution allows us to resolve 
smaller radii is important since the baryons often probe just 
the central few percent of a dark matter structure - the lat- 
est observations of galaxy and clusters probe the mass dis- 
tribution within one percent of the virial radius, which un- 
til recently was unresolved by numerical simulations. Forth - 
coming experi ments, such as VERITAS l|Weekes et. al2002l) 
and MAGIC jFlix. Martinez fc Pradall2004) will probe the 
structure of dark matter halos on even smaller scales by 
attempting to detect gamma-rays from dark matter annihi- 
lation within the central hundred parsecs Q-l%Rviriai ) 
of the Galactic halo fe.g. ICalcaneo-Roldan fc Moorell200(i) . 
These scales are still below the resolution limit of todays 
cosmological simulations, the estimates of the dark matter 
densities in these regions are still based on extrapolations 
which introduce large uncertainties 

A simple estimate of the scaling of A'^ with time 
shows remarkable progress over and above that predicted 
by Moore's law. The first computer simulations used of the 
order lO'^ partic les and force resolution of the order of the 
half mass radn (|Peeblejll970l) . Today we can follow up to 
10* particles with a resolution of 10~^ of the final struc- 
ture. The increase in resolution is significantly faster than 
predicted by Moore's law since equally impressive gains in 
performance have been due to advances in software. 

We are finally at the stage whereaby the dark matter 
clustering is understood at a level where the uncertainties 
are dominated by the influence of the baryonic component. 
It is therefore a good time to review and compare existing 
results from different groups together with a set of new simu- 
lations that we have carried out that are the state of the art 
in this subject and represent what is achievable with sev- 
eral months of dedicated supercomputer time. For certain 
problems, such as predicting the annihilation flux discussed 
earlier, it would be necessary to significantly increase the 
resolution. This is not possible with existing resources and 
new techniques should be explored. We begin by presenting 
our new simulations in Section 2. Section 3 discusses con- 
vergence tests and the asymptotic best fit density profiles. 
In Section 4 we compare our results with recently published 
results from four other groups mentioned above. 



2 NUMERICAL EXPERIMENTS 

Table gives an overview of the simulations we present in 
this paper. With up to 25 x 10^ particles inside the virial 
radius of a cluster and an effective 10® timesteps, they are 
among the highest resolution ACDM simulations performed 
so far. They represent a major investment of computing 
time, the largest run was completed in about 10® CPU hours 
on the zBox supercomputer ^ . 



^ http:/ /www-theorie.physik.unizh.ch/~stadel/zBox/ 



2.1 N-body code and numerical parameters 

The simulations have been performed using a new version of 
PKDGRAV, written by Joachim Stadel and Thomas Quinn 
JStade<l200lh . The code was optimised to reduce the com- 
putational cost of the very high resolution runs we present 
in this paper. We tested the new version of the cod e by re- 
runn ing the "Virgo cluster" initial conditions 4Moore et alJ 
Il998h . We confirmed that density profile, shape of the clus- 
ter and the amount of substructure it contains is identi- 
cal to that obtaine d with the original code presented in 
IChignaet all dlOgsTl . 

Individual time steps are chosen for each particle pro- 
portional to the square root of the softening length over the 
acceleration, /Sti — r\yj t/ai. We use r] = 0.2 for most runs, 
only in run D'dlt we used larger timesteps 77 = 0.3 for com- 
parison. The node-opening angle is set to 6^ = 0.55 initially, 
and after z = 2to 9 — 0.7. This allows higher force accuracy 
when the mass distribution is nearly smooth and the rela- 
tive force errors can be large in the treecode. Cell moments 
are expanded to fourth order in PKDGRAV, other treecodes 
typically use just second or first order expansion. The code 
uses a spline softening length e, forces are completely New- 
tonian at 2e. In Table 0eo is the softening length at 2; = 0, 
^max is the maximal softening in comoving coordinates. In 
most runs the softening is constant in physical coordinates 
from z — 9 to the present and is constant in comoving coor- 
dinates before, i.e. 6max — lOeo ■ In runs C9 and FQctti the 
softening is constant in comoving coordinates for the entire 
run, in run F9ft the softening has a constant physical length 
for the entire run. 



2.2 Initial conditions and cosmological parameters 

We adopt a ACDM cosmological model with parameters 
from the first year WMAP r esults: Qa = 0.7 32, Qm = 
0.268, as = 0.9, h = 0.71, SSpergel et al.ll2003l) . The ini- 
tial conditions are generated with the GRAFIC2 package 
jBertschingeJbOOll) . The starting redshifts Zi are set to the 
time when the standard deviation of the density fiuctuations 
in the refined region reaches 0.2. 

First we run a parent simulation; a 300^ particle cu- 
bic grid with a comoving cube size of 300 Mpc (particle 
mass rup — 3.7 x IO'^'^Mq, force resolution eo = lOOfcpc, 
Emax = I Mpe). Then we u se the friends-of-friends (FoF) 
algorithm jDavis et al.l ll985') with a linking length of 0.164 
mean interparticle separations to identify clusters. We found 
39 objects with virial masses above 2.3 x 10^'' M©. We se- 
lected six of these clusters for resimulation, discarding ob- 
jects close to the periodic boundaries and objects that show 
clear signs of recent major mergers at 2; = 0. We label the six 
cluster with letters AXo F according to their mass. It turned 
out that two of the clusters selected in this way (runs A and 
C) have ongoing major mergers at 2; = (i.e. two clearly 
distinguishable central cores) , which is not evident from the 
parent simulation due to lack of resolution. These clusters 
were evolved slightly into the future to obtain a sample of 
six 'relaxed' clusters. 

For re-simulation we mark and trace back the particles 
within a cluster's virial radius to the initial conditions. All 
particles which lie within a 4 Mpc (comoving) thick region 
surrounding the marked particles in the initial conditions are 
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Table 1. Parameters of simulated cluster halos 
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also added to the refinement region. Tiiis ensures tiiat tiiere 
is no pollution of heavier particles within the virial radius of 
the resimulated cluster. Typically one third or one quarter 
of the refinement particles ends up within the virial radius. 
To reduce the mass differences at the border of the refine- 
ment region we define a 5 Mpc thick 'buffer region' around 
the high resolution region, there an intermediate refinement 
factor of 3 or 4 in length is used. The final refinement fac- 
tors are 6, 9 and 12 in length, i.e., 216, 729 and 1728 in 
mass, so that the mass resolution is rup = 2.14 X IO^Mq m 
the highest resolution run. We label each run with a letter 
indicating the object and number that gives the refinement 
factor in length. To reduce the mass differences at the bor- 
der of the refinement region we define a 5 Mpc thick 'buffer 
region' around the high resolution region, there an interme- 
diate refinement factor of 3 or 4 in length is used. 



3 ACDM CLUSTER PROFILES 
3.1 Profile Convergence Tests 

Numerical convergence tests show that with sufficient 
timesteps, force accuracy and force resolution the radius a 
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2001^ 
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compare different mass resolution simulations of the same 
object to determine the resolved r adius. The r e sultin g 
radii scale with A'^"'' *^ accordi ng to | Pm ver_gt aP ||2Qq3)) 
igavashi ct al. (2003) and Fukush ig^ et al.l (l20oJrbiit only 
with j V~^^ ^ in the tests inlMoore et al.l Jl998t)lGhigna etaO 
\200(t and lR,eed etH (|2£ 



2.3 Measuring density profiles 

We define the virial radius r^ir such that the mean den- 
sity wi thin rvir is nSQlf^pcrit = 98.4pcrit for the adopted 
model llEke. Cole fc Frenklll99i ). We use 30 spherical bins 
of equal logarithmic width, centered on the densest region of 
each cluster using TIPSY ^ We confirmed that using triaxial 
bins adapted to the shape of the isodensity surfaces (at some 
given radius, we tried 0.1, 0.5 and Irvir) does n ot change the 
form o f the density profile, in agreement with l.Ting fc Sutol 
i2002 Y For simplicity and easier comparison to other results 
we will present only profiles obtained using spherical bins. 
Data points are plotted at the arithmetic mean of the cor- 
responding bin boundaries; the first bin ends at 1.5 kpc, the 
last bin at the virial radius. 



3.1.1 Mass resolution 

The finite mass resolution of N body simulations always 
leads to two body relaxation effects, i.e. heat is transported 
into the cold halo cores and they expand. It is not obvious 
that better mass resolution reduces the effects of two body 
relaxation, since in hierarchical models the first resolved ob- 
jects always contain just a few particles and with higher res- 
olution these first objects form earlier, i .e. they are denser 
and more affected by r elaxation effects ijMoore et alJl200ll : 
iBinnev fc Knebell2002h . Estimates of relaxation based on 
following the local phase-space density in simulations show 
that the amount of relaxation can be reduced with better 
mass resolution, but the average degree of relaxation scales 
roughly like N'"'^ much slower than expected from th e 
relaxation time of the final structure (jPiemand et alJl2004^ . 
This confirms the validity of performing convergence tests 
in A'^, but one has to bear in mind that convergence can be 
quite slow. 

We checked a series of resimulations of the same cluster 
(D) for convergence in circular velocity, mass enclosed and 



2 TIPSY is available form the Uni- 

versity of Washington N-body group: 

|http://www-hpcc. astro, washington.edu/tools/tipsy /tipsy.html| 



^ Convergence within 10% in cumulative mass is the same as 
convergence in circular velocity with a tolerance of 5% 
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Figure 1. Numerical convergence tests for tlie cluster profiles: Panel (a): Density profiles of cluster D resolved with A'^vir = 205fc, 1.8M, 6M 
and 14M particles. Panel (b): Logarithmic slope for the profiles from (a). Panel (c); Density profiles of cluster F simulated w ith different 
numerical parameters: F9ft used 4096 fixed timesteps and constant e in physical coordinates as in lFukushige et alj i2004l) . F9cm and 
_F9 used adaptive timesteps 0.2-^/ €{z)/a with comoving softening in F9 and mixed comoving/physical softening in _F9 (emax = lOeo). 
Panel (d): Logarithmic slope for the profiles from (c). 



density. Outside of the converged radii the values must be 
within 10% of the reference run D12. Table |5| shows the 
measured converged radii. 

(i) Convergence is slow, roughly oc N~^^^. Therefore a 
high resolution reference run should have at least 8 times 
as many particles. Between run D9 and D12 the factor is 
only 2.37. Using D12 to determine the converged radii of D9 
gives radii that are ab out a factor two too small (Table I^J. 
iFukushige et al.1 (|2Qq3) compare runs with A'^^,ir = 14 x lO*" 
and A'vir = 29 x 10®. At radii where both runs have similar 
densities it is still not clear if the simulations have converged, 
even higher resolution studies are needed to demonstrate 
this. 

(ii) If one sets the force resolution to one half of expected 
resolved radius, then it is not surprising to measure a re- 



solved radius close to the expected value. With this method 
one can demonstrate almost arbitrary convergence criteria, 
as long as they overestimate rconv Therefore convergence 
tests in A'^ should be performed with small softenings (high 
force resolution). Runs D3h, D6h and D12 all have eo = 1.8 
kpc, their converged radii scale like the mean interparticle 
separation N~^^^. I n run D6 ep = 3.2 kpc is close to the 
'optimal value' from lPower et all lf2003) , and the converged 
radii are larger than in D6h (see Figure 

(iii) Different small scale noise in the initial conditions 
leads to different formation histories. Therefore the shape 
and the density profile can differ even at radii were all runs 
have converged. For example between r=10 kpc and 320 kpc 
the densities in run D9 are about 7% higher than in run D12. 
Therefore the densities in D9 are within 10% from those of 
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Figure 2. Ratios of tlic mass enclosed in low resolution runs to 
mass enclosed in tlie higli resolution run D12. By comparing runs 
with equal softening (smaller than one third of the convergence 
scale) like D3h and D6h one finds that the resolved radii scale 
like r oc N~^^^. A larger softening (see run D6) can increase the 
converged scales and change this scaling. 



Table 2. Convergence radii measured by comparing with run 
D12. The numbers in the run labels are oc N^i'^, at fixed force 
resolution we get r oc N~^/^ (bold values). Question marks indi- 
cate that a run with much better mass resolution than D12 would 
be needed to measure this convergence radii reliably. Stars indi- 
cate estimated radii assuming a convergence rate of r oc N~^/^. 
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D12 quite early. If one rescales p in this range rio%p of D9 
grows from 2.2 kpc to 4.6 kpc. 

Extrapolating rconv oc N~^^^ to our highest resolution 
runs gives the values on the last two lines of Table |21 Note 
that this is just an extrapolation, it is not clear that this 
scaling is valid down to this level, only larger simulations 
could verify this. To be conservative we assume the limit 
due to mass resolution to be 9 kpc for the '9-series' of runs, 
and 6.8 kpc for run D12. The for ce resolution sets another 
limit at about 3eo llMoore et all . 1998 Ghigna ct al. 2000). 
We give the larger of the two limits as the trusted radius in 
Table [U 
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Figure 3. The triangles show the timestep criterium r]y'e{z)/a 
as a function of radius for run D9 at z = 0. The dashed line is for 
run D9lt, which has rj = 0.3, and the long dashed line for run D12. 
The o pen squares give 15(At/to)^/^tcirc(''uir) form lPower et all 
l200.'f> . the circles are the circular orbit timescale 2iTr /vdrci''')- 
Lines without symbols show = l/(-\/ Gp{< r)15). The 

two horizontal lines are the timesteps and 15(At/to)^''^tcirc(''i>ir) 
for run F9ft. 

3.1.2 Force and time resolution 

Finite timesteps and force resolution also sets a limiting ra- 
dius/density that a run can resolve. We use multistepping, 
individual timesteps for the particles that are obtained by 
dividing the main timestep (usually to/200) by two until 
it is smaller than r)^ e{z)/a, where a is the local accelera- 
tion. Our standard choice is rj = 0.2 and e{z = 0) between 
O.OOlrvir and 0.0022rvir, which is chosen to be less than one 
third of the resolution limit expected from the finite mass 
resolution, e is constant in physical length units since 2 = 9 
and comoving before that epoch. Here we argue that the 
resolution limit imposed by this choice of multistepping lie 
well below the scale affected by finite mass resolution. 

In run D9lt the number of timesteps was reduced by 
using rj = 0.3, at equal force resolution as in D9. Run F9cm 
had a constant comoving softening during the entire simula- 
tion, in run F9ft the softening is physical and the timesteps 
are fixed At = to/4096 and equa l for all particles (i.e. t he 
same numerical parameters as in iFukushige et al.l ll2004l) ). 
The density profiles are very similar (Figure Panel c), 
there is no significant difference above the mass resolution 
scale of 9 kpc. There is a small difference in the inner pro- 
file of F9 compared to -F9/t and F9cm, at large z this run 
has larger e and therefore larger timesteps than F9cm. So 
it is possible that runs with our standard parameters have 
slightly shallower density profiles at the resolution limit than 
runs with entirely comoving softening, or runs with a suffi- 
ciently large number of fixed timesteps. However run F9cm 
takes twice as much CPU time as run F9 and run F9/t 
three times more, therefore we accept this compromise. 
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Figure 121 shows the timestep criterion r)^ e{z)/a as a 
function of radius at z = for runs D9 (triangles, solid line), 
D%lt (dashed) D12 (long dashed) and for F'dft (horizontal 
line). Particles near the cluster centre must take timesteps 
below 2 X lQ~'^to, i.e. their timesteps are to/200 x — 
to/6400. According to IPower et al.l ll2003l) the resolution 
limit due to finite timesteps tts is where the circular velocity 
(circles) equals 15(At/to)^^^tcirc(f„ir) (open squares). This 
radius is indeed close to that where the circular velocities 
and densities start to differ, however for run DQlt this esti- 
mate is even a bit too conservative, since the density (and 
also Wcirc) profiles of DQlt and -D9 agree down to at least 
O.OOSrvir. This suggest that about 15 timesteps per local 
dynamical time are sufficient for the simulations presented 
here. Note that this is probably not a general condition for 
all cosmological simulations: Other codes seem to require 
different convergence co nditions than those we present in 
this paper. For example. iFukushige et al.l (1200411 found that 
their runs converge down to O.OOSrvir even with only 2048 
fixed timesteps, which corresponds to only eight timesteps 
per dynamical time at this radius. 



3.2 Density Profiles 

In this section we present the profiles of the six high reso- 
lution runs: A9,B9,C9,D12,£;9,F9cm. The output at z = 
was used, except for clusters A9 and C9 which had a re- 
cent major merger * and the core of the infalling cluster is 
at about 0.02rvir in A9 and at O.lrvir in C9. These cores 
spiral in due to dynamical friction and in the 'near' future 
both clusters have a regular, 'relaxed' central region again. 
Therefore we use outputs ai z = —0.137 (-1-2.1 Gyr) for run 
A9 and z = -0.167 (+2.6 Gyr) for C9. 



3.3 Two parameter fits 

Figure |1] shows the density profiles of the six different 
clusters. We also show best fits to functions previously 
proposed in the literature that have asy mptotic central 
slopes of -1 i Navarro. Frenk fc Whitelll996l :NFW) and -1.5 
jMoore et alJ 1999 ;M99). The fits are carried out over the 



tral slope and propose a profile which curves smoothly over 
to a constant density at very small radii: 



resolved region by minimising the mean square of the rel- 
ative density differences. These two profiles have two free 
parameters, namely the scale radius and the density at 
this radius ps — p(rs). The scale radii of these best fits 
give the concentrations c = Vvii/ra listed in Table |3] The 
residuals are plotted in the top and bottom panels of Fig- 
ure |1] and the rms of the residual are given in Table 13 as 
Anfw and Am99. The residuals are quite large and show 
that neither profile is a good fit to all the simulations which 
lie somewhere in between these two extremes. 



3.4 Three parameter fits 

iNavarro et all J2003l) argue the large residuals of NFW and 
M99 fits are evidence against any constant asymptotic cen- 



ln{pN{r)/ps) = (-2/q:n) [{r/rs 



(1) 



This function gives a much better fit to the simulations, 
see the dashed dotted lines in Figure but this should be 
expected since there is an additional third free parameter 
On, while the NFW and M99 profiles only have two free 
parameters, qn determines how fast the profile ^ turns 
away from an power law near the centre. INavarro et alJ 
l|2003l) found that q:n is independent of halo mass and 
Qn ~ 0.172 ±0.032 for all their simulations, including galax- 
ies and dwarfs. The mean and scatter of our six high reso- 
lution clusters is q:n ~ 0.186 ±0.037. (Excluding cluster C9 

yields qn = 0.174 ± 0.025). 

We also show fits to a general a/37-profile JZhadll99d) 
{pa, subscript 'G' stands for 'general') that asymptotes to 
a central cusp p{r) cx r~^: 



pair) 



(r/r,)T(l±(r/r,)")W-^)/° 



(2) 



We fix the outer slope (3 — 3 and the turnover parameter a = 
1. For comparison the NFW profile has (a,/?, 7) = (1,3, 1), 
the M99 profile has (a,/3,7) = (1.5, 3, 1.5). We fit the three 
parameters 7, Vs and ps to the data and find that this cuspy 
profile also provides a very good fit to the data. The best fit 
values and rms residual are listed in Table 3 and we find a 
mean slope of 7 = 1.16 ± 0.14. 

Using a sharper turnover a = 1.5 makes the fits slightly 
worse (the average of Ac is about 20 percent larger) and 
the best fit inner slopes are somewhat steeper 7 = 1.31 ± 
0.11. We also made some attempts with fitting procedures 
wher e a ot or bot h a and /3 are also free parameters. Like 
Klvp in et alJ (120011) we found strong degeneracies, i.e. very 
different combinations of parameter values can fit a typical 
density profile equally well. Therefore we only present results 
from the fits with fixed a and /3 parameters in this paper. 

The fitting functions Q and ^ fit the measured den- 
sity profiles very well over the whole resolved range. Func- 
tion 10 is even a relatively good approximation below the 
resolved scale: For example if one is extremely optimistic 
about rrcsoivcd in run D6 and uses rrcsoived = 2.8 kpc in- 
stead of 13.5 kpc one gets qn = 0.0203, cn = 7.1 and 
An ~ 0.127, while the generalised fit is now clearly worse: 
7G = 0.99, CG = 3.6 and Ac = 0.216. Also note that the 
residuals near rrcsoivcd are very small or positive for (0, i.e. 
the measured density is as large as the fitted value. But at 
rresoived it is possiblc that the measured density is slightly 
too low since in this region the numerical limitations start 
to play a role. If extrapolation beyond the converged radius 
is necessary it is not clear which profile is a safer choice. 
We agree with Navarro ct al. (200j) that all simple fitting 
formula have their drawbacks, that direct comparison with 
simulations should be attempted whenever possible and that 
much higher resolution simulations are needed to establish 
(or exclude) that CDM halos have divergent inner density 
cusps (as predicted in iBinnev...2003i'l . 



^ An mpeg movie of the formation of 
cluster C9 can be downloaded from 

|http://www-theorle.physlk.unizh.ch/~dlemand/cIusters/| 



3.5 Maximum inner slope 

The results from the last section suggest that profiles with 
a central cusp in the range 7 = 1.16 ± 0.14 provide a 
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Figure 4. Density profiles of the six clusters in our sample, clusters B to F are shifted downwards for clarity. Clusters are ordered by 
mass form top to bottom. Profiles of cluster A and C are shown at redshifts —0.14 and —0.17, i.e. when they have reached a 'relaxed' 
state with one well defined centre. Best fit NFW and M99 profiles and residual are shown, obtained by minimising the squares of the 
relative density differences. 

Table 3. Density profile parameters. A is the root mean square of (p — Pflt)/p for the four fitting functions 
used. 



Run 


cnfw 


Anfw 


CM99 


Am99 


7G 


CG 


Ag 


ON 


CN 


An 


A9 


5.7 


0.10 


1.7 


0.21 


1.16 


3.9 


0.057 


0.167 


4.2 


0.033 


B9 


4.2 


0.16 


1.5 


0.13 


1.29 


2.1 


0.083 


0.141 


2.6 


0.093 


C9 


7.6 


0.09 


3.0 


0.26 


0.92 


8.7 


0.081 


0.247 


7.2 


0.068 


D3h 


7.4 


0.17 


3.9 


0.13 


1.42 


4.0 


0.103 


0.175 


7.3 


0.101 




7.9 


0.11 


3.8 


0.13 


1.17 


4.6 


0.089 


0.206 


7.2 


0.081 


De 


7.9 


0.12 


3.8 


0.16 


1.25 


5.4 


0.101 


0.193 


7.2 


0.097 


D9 


8.8 


0.12 


3.9 


0.12 


1.21 


6.2 


0.096 


0.190 


7.8 


0.087 


Dm 


8.7 


0.12 


3.8 


0.12 


1.20 


6.2 


0.098 


0.191 


7.7 


0.087 


D12 


8.4 


0.12 


3.1 


0.14 


1.25 


4.5 


0.066 


0.174 


6.9 


0.051 


E9 


7.4 


0.12 


3.0 


0.10 


1.25 


4.5 


0.072 


0.176 


6.2 


0.069 


F9 


6.9 


0.06 


3.0 


0.14 


1.02 


6.7 


0.054 


0.224 


6.5 


0.048 


F9cm 


7.3 


0.06 


3.1 


0.14 


1.10 


6.2 


0.055 


0.212 


6.6 


0.057 


F9ft 


7.2 


0.05 


3.1 


0.16 


1.05 


6.6 


0.043 


0.218 


6.5 


0.045 
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Figure 5. Same as Figure E| but with fitt ing functions tliat liave one additional free parameter. The dashed dotted lines show the 
profile proposed by I Nava rro et al.lEoO.'j) . The dashed lines show a general a/37-profile . We fitted the inner slope 7 to the data 
and used fixed values for the outer slope /3 = 3 and turning parameter a = 1. 7 = 1 corresponds to the NFW profile. The fit parameters 
and rms of the residuals are given in Table El 



good approximation to the inner dens i ty pro files of ACDM 
halos. But Figure 4 in iNavarro et al.l Il2003() seems to ex- 
clude our mean value for more than half of their cluster 
profiles. This is not totally inconsistent, but a hint for a 
mild discrepancy that we will try to explain: In principle 
the mass inside the converged radius limits the inner slope: 
7max ~ 3(1 — p{r)/ p(,r)). This is true if both the density 
and the cumulative density are correct down to the resolved 
scale. But up to now the central density of a simulated pro- 
file always increased with better numerical resolution, so it 
is likely that also todays highest resolution simulations un- 
derestimate the dark matter density near the centre. This 
means that cumulative quantities like i;circ(r), M(< r) and 
p(< r) tend to be too low even at radii wh ere the density 
has co nverged. The converged radii used in iNavarro et alJ 
i2003l) are close to the radius where the circular velocity is 
within 10 percent of a higher resolution run, while the den- 
sity converges further in at about O.Grconv jHavashi ct al. 
1200311 . If we assume that this is also true for their highest 
resolution runs then p(< r) oc Vcirc{r)'^ is up to 20 percent 



too low, while the error in p(r) is much smaller. This raises 
the values for 7max by about 0.2 ~ 0.3 and our mean value 
7 — 1.16 is not excluded by any of their clusters anymore. 
If the convergence with mass resolution is not as fast as 
J'conv oc A'^"^'*^ but rather rconv oc A'^^^''^, see Section l3.ll 
then the maximum inner slopes could have even larger er- 
rors. 



4 COMPARISON WITH OTHER GROUPS 

Recently, several groups have published simulations of dark 
matter clusters in the concordance cosmological model. 
These authors kindly suppli ed their density p r ofiles and we 
show the comparison here. iFukushige et al.1 J2004 CF03'') 
simulated four ACDM clusters with 7 to 26 million particles 
using a Treecode and the GRAPE hardware. These authors 
also used the GRAFIC2 softwa re tBcrtschingc r 2001) to gen- 
erat e their initial condit ions. iHavashi et al.l ll2003i) ('H03'') 
and INavarro et alJ (1200311 presented eight clusters resolved 
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with up to 1.6 milli on particles within r 2oo simulated with 
the GADGET code JSpringel et al.ll200j) . the method used 
to gen e rate the initial condition s is described in lPower et alJ 
J2003h . iTasitsiomi et all (|20o3)('T03') simulated six clus- 
ters with up to 0.8 million particles wi thin riso using the 
adaptive refinement tree code ART (iKravtsov et alJll997l) 
and a technique for setting up mult i- mass initial condi- 
tions described in iKlvpin et al.l J200ll) . IWambsganss et alJ 
J20^)('W03') present a cosmological simulation without 
resimulation of refined regions, i.e. constant mass resolution 
(1024'^ particles in a 320 /i^^Mpc box). The four most mas- 
sive clusters in this cube are resolved within 0.5 to 0.9 mil- 
fion particles. This simulation was performed with a Tree- 
Particle-Mesh (TPM) code (Bode fc QstrikedboOSh with a 
softening of 3.2 h~^kpc. 

In Figure|S]we show these data along with the new sim- 
ulations presented in this paper. We plot the density profiles 
and the logarithmic slopes of the clusters all normalized at 
the radius such that the circular velocity curve peaks rvcmax 
and to p(< rvcmax). This corresponds to the radius at which 
dlogp/dlogr = —2. We plot the curves to the "believable" 
radius stated by each group and down to about 0.01 r-vir for 
W03. 

The density profiles are reassuringly similar. Further- 
more, the scatter is small, the standard deviation of all pro- 
files is roughly ±0.15 in the logarithmic gradient at small 
radii (0.01 — 0.5rvcmax). Table 0] lists the measured slopes 
at different radii. There is no valu e at 3%rvcmax for the 
cluste r from lTasitsiomi et al.l i200'j) and lWambsganss et alJ 
i2004l) because this is below their quoted resolution limit. 
Most values agr ee within the scatter, the profiles from 
iTasitsiomi et all l|^04) are steeper when compared at 0.01 
and 0.03 rvir = ?'98.4, but within the scatter at 3%rvcmax. 
This could be due to different halo selection. The majority 
of their clusters are not isolated but in close pairs or triplets. 
In a close pair the density falls slower with radius to 98.4 
pcrit, SO rvir = '"98.4 is further out as in a isolated cluster 
with similar inner profile. Among the samples of isolated 
clusters (F03; H03; W03 and our clusters) there is a small 
trend at 0.01 Tvir towards steeper slopes with better mass 
resolution. This could indicate that some numerical flatten- 
ing of the profiles is still present at 0.01 rvir in the lower 
resolution clusters. 



5 SUMMARY 

We have carried out a series of six very high resolution cal- 
culations of the structure of cluster mass objects in a hi- 
erarchical universe. The clusters contain up to 25 million 
particles and have force softening as small as 0.1%rvir. 

A convergence analysis demonstrates that for our 
Treecode with our integration scheme, the radius beyond 
which we can trust the density profiles scale according to 
the mean interparticle separation. In the best case we reach 
a resolution of about 0.3%rvir. 

Neither of the two parameter functions, the NFW and 
M99 profiles, are very good fits over the whole resolved range 
in most clusters. One additional free parameter is needed 
to fit all six c l usters : The asymptotically flat profile from 
iNavarro et al.l i2003h and an NFW profile with variable in- 
ner slope provide much improved fits. The best fit inner 



Table 4. Logarithmic slopes (absolute values) of our six high 
resolution cluster density profiles. Line (a) gives the averages and 
scatter, (b)-(c) are average slopes from other groups (see text for 
details). 



l%rvir 3%rvir 3%rv cmax 9%rv CI 



A9 1.22 1.36 1.24 1.64 

B9 1.33 1.43 1.21 1.63 

C9 1.24 1.21 1.25 1.26 

D12 1.28 1.54 1.32 1.58 

E9 1.31 1.44 1.41 1.62 

F9cm 1.19 1.47 1.22 1.43 



a) A-F 1.26±0.05 1.41±0.11 1.28±0.08 1.53±0.15 

b) F03 1.25±0.05 1.52±0.06 1.33±0.15 1.54±0.15 

c) H03 1.18±0.13 1.38±0.14 1.23±0.17 1.50±0.14 

d) T03 1.50±0.14 1.79±0.07 - 1.56±0.12 

e) W03 1.11±0.04 1.41±0.13 - 1.35±0.06 



avg.(a-o) 1.26 1.50 - 1.49 

avg.(a-c) 1.23 1.44 1.28 1.52 



slopes are 7 = 1.16 ± 0.14. Below the resolved radius the 
two fitting formulas used are very different. Future simula- 
tions with much higher resolution will show which one (if 
either) of the two is still a good approximation on scales of 
0.1%r-vir and smaller. 

We compare our results with simulations from other 
groups who used independent codes and initial conditions. 
We find a good agreement between the cluster density 
profiles calculated with different algorithms. From 0.03 — 
0.5rvcmax the scatter in the profiles is nearly constant and 
equal to about 0.17 in logarithmic slope. At one percent of 
the virial radius (defined such that the mean density within 
rvir is 178n5/^ Pcrit = 98.4pc:rit) the slope of the density pro- 
files is 1.26 ±0.16. 
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